PALEVO 1631-0683 Elsevier S1631-0683(19)30186-1 10.1016/j.crpv.2019.10.008 Research article Palaeontology, Systematics and Evolution (Evolutionary Patterns, Macroevolution, Neontology) Early tetrapodomorph biogeography: Controlling for fossil record bias in macroevolutionary analyses Biogéographie des premiers tétrapodomorphes : contrôle du biais d’enregistrement des fossiles dans les analyses de macroévolution Folie Annelise Buffetaut, Eric Bardet Nathalie Houssaye Alexandra Gheerbrant Emmanuel Laurin Michel Gardner Jacob D. a 1 Surya Kevin b 1 Organ Chris L. organ@montana.edu a Department of Earth Sciences, Montana State University, 59717 Bozeman, MT, USA Department of Earth Sciences, Montana State University Bozeman MT 59717 USA Honors College, Montana State University, Bozeman, 59717 MT, USA Honors College, Montana State University Bozeman MT 59717 USA

Contributed equally to this work.

 18 7 S1631-0683(19)X0008-1 Palaeobiology and palaeobiogeography of amphibians and reptiles: An homage to Jean-Claude RagePaléobiologie et paléobiogéographie des amphibiens et reptiles : un hommage à Jean-Claude Rage 699 709 © 2019 Académie des sciences. Published by Elsevier B.V. All rights reserved. 2019 Académie des sciences Full (PDF)

The fossil record provides direct empirical data for understanding macroevolutionary patterns and processes. Inherent biases in the fossil record are well known to confound analyses of this data. Sampling bias proxies have been used as covariates in regression models to test for such biases. Proxies, such as formation count, are associated with paleobiodiversity, but are insufficient for explaining species dispersal owing to a lack of geographic context. Here, we develop a sampling bias proxy that incorporates geographic information and test it with a case study on early tetrapodomorph biogeography. We use recently-developed Bayesian phylogeographic models and a new supertree of early tetrapodomorphs to estimate dispersal rates and ancestral habitat locations. We find strong evidence that geographic sampling bias explains supposed radiations in dispersal rate (potential adaptive radiations). Our study highlights the necessity of accounting for geographic sampling bias in macroevolutionary and phylogenetic analyses and provides an approach to test for its effect.

Le registre fossile fournit des données empiriques directes pour la compréhension des patrons et processus macro-évolutifs. Il est bien connu que des biais inhérents à l’enregistrement fossile peuvent brouiller les analyses de ces données. Des variables de substitution des biais d’échantillonnage ont été utilisées comme covariables dans des modèles de régression pour vérifier la présence de tels biais. Des approximations comme le dénombrement des formations sont associées à la paléobiodiversité, mais elles sont insuffisantes pour expliquer la dispersion des espèces en raison d’un défaut de contexte géographique. Ici, nous élaborons un substitut de biais d’échantillonnage qui incorpore de l’information géographique, et nous le testons sur une étude de cas de biogéographie sur les premiers tétrapodomorphes. On utilise des modèles phylogéographiques bayésiens récemment développés et un nouvel arbre de synthèse des premiers tétrapodomorphes pour estimer les taux de dispersion et les localisations d’habitats ancestraux. Nous obtenons des preuves solides de ce que le biais d’échantillonnage géographique explique les radiations supposées dans le taux de dispersion (radiations adaptatives potentielles). Notre étude met l’accent sur la nécessité de tenir compte des biais de l’échantillonnage géographique dans les analyses macro-évolutives et phylogénétiques, et fournit une approche pour tester leur effet.

 Sampling bias, Fossil record, Biogeography, Phylogenetics, Macroevolution, Tetrapodomorph water–land transition Biais d’échantillonnage, Registre fossile, Biogéographie, Phylogénétique, Macro-évolution, Sortie de eaux chez les Tétrapodomorphes presented Handled by Michel Laurin Introduction

Our understanding of macroevolutionary patterns and processes are fundamentally based on fossils. The most direct evidence for taxonomic origination and extinction rates come from the rock record, as do evidence for novelty and climate change unseen in data sets gleaned from extant sources. There are no perfect data sets in science; there are inherent limitations and biases in the rock record that must be addressed when we form and test paleobiological hypotheses. For instance, observed stratigraphic ranges of fossils can mislead inferences about diversification and extinction rates (Raup and Boyajian, 1988 and Signor and Lipps, 1982). Observed species diversity is also known to increase with time due to the preferential preservation and recovery of fossils in younger geological strata—referred to as “the Pull of the Recent” (Jablonski et al., 2003). Large and long-surviving clades with high rates of early diversification tend to result in an illusionary rate slow-down as diversification rates revert back to a mean value—referred to as “the Push of the Past” (Budd and Mann, 2018). Paleobiologists test and account for these biases when analyzing diversification and extinction at local and global scales (Alroy et al., 2001, Benson and Butler, 2011, Benson and Upchurch, 2013, Benson et al., 2010 and Benton et al., 2013; Foote, 2003, Jablonski et al., 2003, Koch, 1978, Lloyd, 2012, Sakamoto et al., 2016a and Sakamoto et al., 2016b). These bias-detection and correction techniques include fossil occurrence subsampling (Alroy et al., 2001, Close et al., 2019, Dunne et al., 2018, Jablonski et al., 2003 and Lloyd, 2012); correcting origination, extinction, and sampling rates using evolutionary predictive models (Foote, 2003); the use of residuals from diversity-sampling models (Benson and Upchurch, 2013, Benson et al., 2010 and Sakamoto et al., 2016b); and the incorporation of sampling bias proxies as covariates in regression models (Benson and Butler, 2011, Benson et al., 2010, Benton et al., 2013 and Sakamoto et al., 2016a). Benton et al. (2013), studying sampling bias proxies, demonstrated that diversity through time closely tracks formation count (Benton et al., 2013). However, case studies in England and Wales suggest that proxies for terrestrial sedimentary rock volume (such as formation count) do not accurately explain paleobiodiversity, particularly if the fossil record is patchy (Dunhill et al., 2013, Dunhill et al., 2014a and Dunhill et al., 2014b). Marine outcrop area and paleoecological-associated facies changes are, however, associated with shifts in paleobiodiversity (Dunhill et al., 2013 and Dunhill et al., 2014b). Moreover, Benton et al. (2013) argue that the direction of causality between paleobiodiversity and formation count is unclear; there may be a common cause to explain their covariation, such as sea level (Benton et al., 2013). Nonetheless, formation count is a widely-used sampling bias proxy in phylogenetic analyses of macroevolution (O’Donovan et al., 2018, Sakamoto et al., 2016a, Sakamoto et al., 2016b, Tennant et al., 2016a and Tennant et al., 2016b). The advent of computational modeling approaches, particularly phylogenetic comparative methods, has made it easier to include proxies, like formation count, into models. Additional sampling bias proxies used in these studies include occurrence count, valid taxon count, and specimen completeness and preservation scores. Absent from these proxies is geographic context, which could confound many types of macroevolutionary analyses.

Despite advancements made in understanding the origin and evolution of early tetrapodomorphs, biogeographical studies are hindered by the incompleteness of the early tetrapodomorph fossil record. For example, “Romer's Gap” represents a lack of tetrapodomorph fossils from the end-Devonian to mid-Mississippian, a period crucial for understanding early tetrapodomorph diversification. Recent collection efforts recovered tetrapodomorph specimens from “Romer's Gap”, suggesting that a collection and preservation bias explains this gap (Clack et al., 2017 and Marshall et al., 2019). In addition, a trackway site in Poland demonstrates the existence of digit-bearing tetrapodomorphs 10 million years before the earliest elpistostegalian body fossil, showcasing the limitation of body fossils to reveal evolutionary history (Niedźwiedzki et al., 2010). A recent study by Long et al. (2018) leveraged phylogenetic reconstruction of early tetrapodomorphs to frame hypotheses about the origin of major clades, as well as their dispersal patterns, including the hypothesis that stem-tetrapodomorphs dispersed from eastern Gondwana to Euramerica. However, this study did not use phylogenetic comparative methods to estimate ancestral geographic locations or to model dispersal patterns.

Here, we present a phylogeographic analysis of early tetrapodomorphs. Our goals are: (1) to construct a phylogenetic supertree of early tetrapodomorphs that synthesizes previous phylogenetic reconstructions; (2) to estimate the paleogeographic locations of major early tetrapodomorph clades using recently-developed phylogeographic models that account for the curvature of the Earth; and (3) to test for the influence of geographic sampling bias on dispersal rates. Our results indicate that geographic sampling bias substantially confounds analyses of dispersal and paleogeography. We conclude with a discussion about the necessity of controlling for fossil record biases in macroevolutionary analyses.

 Materials and methods Nomenclature

Tetrapoda has been informally defined historically to include all terrestrial vertebrates with limbs and digits (Laurin, 1998). Gauthier et al. (1989) first articulated a phylogenetic definition of Tetrapoda as the clade including the last common ancestor of amniotes and lissamphibians. This definition excludes stem-tetrapodomorphs, like Acanthostega and Ichthyostega. Stegocephalia was coined by E.D. Cope in 1868 (Cope, 1868), but was more recently used to describe fossil taxa more closely related to tetrapods than other sarcopterygians. A recent cladistic redefinition of Stegocephalia includes all vertebrates more closely related to temnospondyls than Panderichthys (Laurin, 1998). Here, we use the definitions of Laurin (1998) for a monophyletic Stegocephalia and of Gauthier et al. (1989) for Tetrapoda, which refers specifically to the crown group. We use Tetrapodomorpha to refer to all taxa closer to the tetrapod crown-group than the lungfish crown-group (Ahlberg, 1998). We additionally use Elpistostegalia (= Panderichthyida) to refer to the common ancestor of all stegocephalians and Panderichthys as well as Eotetrapodiformes to refer to the common ancestor of all tristichopterids, elpistostegalians, and tetrapods (Coates and Friedman, 2010).

 Supertree

We inferred a supertree of 69 early tetrapodomorph taxa from five edited, published morphological data matrices, focusing on tetrapodomorphs whose previously inferred phylogenetic position bracket the water–land transition (Clack et al., 2017, Friedman et al., 2007, Pardo et al., 2017, Swartz, 2012 and Zhu et al., 2017). Since downstream analyses might be sensitive to unequal sample sizes between taxa pre- and post-water–land transition, we did not include several crownward stem-tetrapodomorphs from the original matrices (Supplementary Material). For each matrix, we generated a posterior distribution of phylogenetic trees using MrBayes 3.2.6 (Ronquist et al., 2012b). In each case, we ran two Markov chain Monte Carlo (MCMC) replicates for 20,000,000 generations with 25% burn-in, each with four chains and a sampling frequency of 1000. We used one partition, except for Clack et al.’s (2017) matrix, which was explicitly divided into cranial and postcranial characters. To time-calibrate the trees, we constrained the root ages and employed a tip-dating approach (Ronquist et al., 2012a). Tip dates (last occurrence) were acquired from the Paleobiology Database (PBDB; https://www.paleobiodb.org/) and the literature (Supplementary Table 2). Root calibrations (minimum and soft maximum age estimates) were collected from the PBDB and Benton et al. (2015). We also used the fossilized birth-death model as the branch length prior (Didier et al., 2017, Didier and Laurin, 2018, Didier et al., 2012, Gavryushkina et al., 2014, Heath et al., 2014, Stadler, 2010 and Zhang et al., 2016). All pairs of MCMC replicates converged as demonstrated by low average standard deviation of split frequencies (< 0.005; Lakner et al., 2008; Supplementary Table 3).

Next, we used the five maximum clade credibility trees (source trees; Supplementary Figs. 1–10) to compute a distance supermatrix using SDM 2.1 (Criscuolo et al., 2006). We then inferred an unweighted neighbor-joining tree (UNJ by Gascuel, 1997) from the distance supermatrix using PhyD* 1.1 (Criscuolo and Gascuel, 2008). The UNJ* algorithm is preferable for matrices based on morphological characters. Unlike most supertree methods, the SDM-PhyD* combination produces a supertree with branch lengths. We rooted the supertree using phytools 0.6.60 (Revell, 2012) by adding an arbitrary branch length of 0.00001 to break the trichotomy at the basal-most node in R 3.5.2 (R Core Team, 2018), designating the dipnomorph Glyptolepis as the outgroup.

We qualitatively compared the supertree topology with the published source trees and Marjanović and Laurin's (2019). We also calculated normalized Robinson-Foulds (nRF) distances (Robinson and Foulds, 1981) using phangorn 2.4.0 (Schliep, 2011) in R to assess the congruency of topologies. In each comparison, polytomies in the supertree or the source tree were resolved in all possible ways using phytools. We then calculated all nRF distances and took an average (Supplementary Table 4). The supplementary materials include a more detailed description of this approach.

 Phylogeography

We obtained paleocoordinate data (paleolatitude and paleolongitude) for 63 early tetrapodomorphs from the PBDB using the GPlates software setting (https://www.gws.gplates.org/). By default, GPlates estimates paleocoordinates from the midpoint of each taxon's age range. Among the 63 taxa sampled, 16 did not have direct paleocoordinate data in the PBDB. For these taxa, we searched for the geological formations and geographic regions within the time range from which they are known and averaged the paleolocations across each valid taxonomic occurrence in the PBDB. If the paleolocation of the formation was not listed in the PBDB, we used published geographic locations of the formations. This level of precision is adequate for world-wide phylogeographic analyses, such as conducted here. Present-day coordinates for these geographic locations were obtained from Google Earth and matched with PBDB entries that date within each taxon's age range (Supplementary Table 5). Four additional taxa, Kenichthys, Koilops, Ossirarus, and Tungsenia, had occurrences in the PBDB but the GPlates software could not estimate their paleocoordinates. For Koilops and Ossirarus, we used all tetrapodomorph occurrences from the Ballagan Formation of Scotland, UK—a formation in which these two taxa are found (Clack et al., 2017). For Kenichthys and Tungsenia, we calculated paleocoordinate data from the GPlates website directly using the present-day coordinates from the PBDB (https://www.gws.gplates.org/#recon-p). This approach did not work for the 16 previously mentioned taxa (Supplementary Table 5). We, therefore, obtained paleocoordinate data from nearby entries in the PBDB that date within each taxon's age range. We excluded the following taxa from our analyses due to the lack of data and comparable entries in the PBDB: Jarvikina, Koharalepis, Spodichthys, and Tinirau. We excluded the outgroup taxon, Glyptolepis, in our analysis to focus on the dispersal trends within early Tetrapodomorpha. We also excluded Eusthenodon and Strepsodus because their high estimated dispersal rates—being reported from multiple continents—masked other rate variation throughout the phylogeny and inhibited our downstream analyses from converging on a stable likelihood. We do, however, discuss their geographic implications in Section 4.

A model that incorporates phylogeny is crucial for paleobiogeographic reconstruction because it accounts for both species relationships and the amount of evolutionary divergence (branch lengths). Using continuous paleocoordinate data, rather than discretely-coded regions, allows dispersal trends to be estimated at finer resolutions. Discretely-coded geographic regions also limit ancestral states to the same regions inhabited by descendant species. However, standard phylogenetic comparative methods for continuous data assume a flat Earth because they do not account for spherically structured coordinates (i.e., the proximity of −179° and 179° longitudes). Recently-developed phylogenetic comparative methods for modeling continuous paleocoordinate data, implemented as the ‘geo’ model in the program BayesTraits V3, overcome this hurdle by “evolving” continuous coordinate data on the surface of a globe (O’Donovan et al., 2018). The model is implemented with a Bayesian reversible jump MCMC algorithm to estimate rates of geographic dispersal and ancestral paleolocations simultaneously. To account for the spheroid shape of the globe, the ‘geo’ model converts latitude and longitude data into three-dimensional coordinates while prohibiting moves that penetrate the inside of the globe. Ancestral states, which are converted back to standard latitude and longitude, are estimated for each node of the phylogeny. The method includes a variable rates model to estimate variation in dispersal rate (Venditti et al., 2011). The ‘geo’ model makes no assumptions about the location of geographic barriers or coastlines, but a study on dinosaur biogeography found 99.2% of mean ancestral state reconstructions to be located within the bounds of landmasses specific to the time at which they occurred (O’Donovan et al., 2018). We ran three replicate independent analyses using the Bayesian phylogenetic ‘geo’ model for 100 million iterations each with a 25% burn-in and sampling every 1000 iterations. We estimated log marginal likelihoods using the Stepping Stone algorithm with 250 stones sampling every 1000 iterations (Xie et al., 2011). We used Bayes factors (BF) to test whether a variable rates model explained the data better than a uniform rate model. Bayes factors greater than two are considered good evidence in support of the model with the greater log marginal likelihood. We compared estimated rate scalars and ancestral states among the three independent variable rates analyses to check for consistency in our results. Rates of dispersal were estimated for each branch by dividing the average rate scalars by the original branch lengths (scaled by time). We assessed the MCMC convergence of all analyses using Tracer 1.7 (Rambaut et al., 2018).

To test for the effect of sampling bias on dispersal rates, we developed a sampling bias proxy that incorporates geographic context: a regional-level formation count. Formation counts are meant to capture multiple biases: uneven global rock exposure, uneven fossil collection and database efforts, and global variation in sediment deposition in environments conducive to preservation. Stage-level (stage-specific) formation count represents the mean number of formations, or distinct rock units, globally known to produce relevant fossils along each terminal branch of a phylogeny. Following the protocol of Sakamoto et al. (2016a) and O’Donovan et al. (2018), stage-level formation counts are calculated by taking the average number of formations known from each geological stage (age) across the globe that encompass the time period between the taxon's tip date and its preceding node. These average stage-level formation counts are weighted by the proportion that each terminal branch length covers each geological stage. For example, if a terminal branch covers two geological stages (e.g., Frasnian and Famennian) at 30% and 70%, respectively, then the formation counts from each geological stage are weighted by those proportions and then divided by the number of geological stages covered: Stage‐level formation count = Frasnian    count × 0 . 3 + Famennian    count × 0.7 2 \[ \text{Stage‐level formation count}=\frac{\text{Frasnian}\text{\hspace{0.28em}count}\times 0.3+\text{Famennian}\text{\hspace{0.28em}count}\times 0.7}{2} \]

Stage-level formation count is not informed by geography; it is a global metric. It is therefore an inadequate proxy if bias has a strong geographic component (e.g., if the majority of formations recorded are from a specific region or if few formations are exposed within a region). The number of fossil-bearing geological formations, accounting for geographic distribution, is expected to be an important confounding bias in the fossil record. We developed a proxy that includes geographic sampling bias. Our approach breaks down stage-level formation count by geographic region. To account for the arrangement of the continents during the Devonian, Carboniferous, and Permian, we recognized five major regions: northern Euramerica (including northeastern Eurasia and central Asia), southern Euramerica (North America, Greenland, and western Europe), western Gondwana (South America and Africa), eastern Gondwana (Antarctica, Australia, and southern Asia), and East Asia (e.g., China). These regions generally resemble traditional bioregionalizations of the Devonian period, but note that regions based on biotic similarities of fossil assemblages are known to change through time (Dowding and Ebach, 2019). Future studies could modify this approach to capture temporal changes in biotic connectivity. For each branch in the phylogeny, we used the average ancestral state and taxon paleolocation estimates to determine if the branch crossed multiple geographic regions. The number of formations within this time window are totaled for every region covered by the branch and then divided by the number of regions covered. For example, if ancestral state estimates at node 1 and 2 are located in eastern Gondwana and southern Euramerica, respectively, then the number of formations recorded in eastern Gondwana, southern Euramerica, and the regions in between (i.e., western Gondwana or northern Euramerica + East Asia) are counted for that geological stage; this total is then divided by the number of geographic regions covered by the entire branch (three for the western Gondwana route and four for the northern Euramerica + East Asia route). If the dispersal path between two consecutive ancestral states does not cross any of the five regions, then the number of formations in the inhabited region is counted alone. Fig. 1 illustrates an example of how this proxy is measured. This results in the average number of formations present along the dispersal path (at geographic region scale) for each branch in the phylogeny. As with stage-level formation counts, the regional-level formation counts are weighted by the proportion that the branch length covers each geological stage. We hypothesize that dispersal rate will inversely correlate with regional-level formation count because we expect that the lack of formations in intermediate regions will lead to inflated dispersal rates. The ‘geo’ model will increase the dispersal rate along a branch to account for the geographic variation observed when there is a lack of intermediate geographic fossil occurrences. This hypothesis can be falsified if high dispersal rates are associated with larger average numbers of formations along dispersal paths. Benton et al. (2013) provide a global sample of tetrapod-bearing rock formations known for each geological stage from the Middle Devonian through the Triassic. We supplemented these lists with stratigraphic units known to produce sarcopterygian fossils entered in the PBDB (collected on December 10th, 2018).

To test for the effect of regional-level formation count bias on dispersal rate, we conducted a non-parametric two-sample, upper-tailed Mann-Whitney U-test using the base package ‘stats’ in R (R Core Team, 2018). This approach ranks all branches of the phylogeny by their regional-level formation count and tests if the branches with lower dispersal rates rank higher on average than branches with higher rates. We define “high” vs “low” dispersal rates based on whether or not they are two standard deviations greater than the average rate across the tree. Due to the vast difference in sample size between the two groups (“high rates”: n = 9, “low rates”: n = 111), we bootstrapped the regional-level formation counts from each group with 100,000 replicates. From this bootstrap analysis, we obtained a 95% confidence interval for the summed ranks of the branches with low dispersal rates (n = 100,000 U-statistic values). The expected U-statistic is 499.5 given the null hypothesis that only 50% of the regional-level formation counts along branches with low rates rank higher than the formation counts with high rates half     of     all     possible     combinations = 9 × 111 2 $ \left(\text{half}\text{\hspace{0.28em}}\text{of}\text{\hspace{0.28em}}\text{all}\text{\hspace{0.28em}}\text{possible}\text{\hspace{0.28em}}\text{combinations}=\frac{9\times 111}{2}\right)$ . A 95% confidence interval of bootstrapped U-statistics that does not include the null expected U-statistic is considered good evidence for higher mean dispersal rates along branches with lower regional-level formation counts. The full dataset and code for the phylogeographic analyses can be requested by email to the corresponding author.

Estimated ancestral states do not identify specific dispersal routes, so we conducted sensitivity analyses to test if the dispersal route chosen for counting formations influenced our results. We conceived of three scenarios for dispersal routes between eastern Gondwana and southern Euramerica or vice versa: (1) a dispersal route through western Gondwana; (2) a route through northern Euramerica and East Asia; and (3) a direct route between eastern Gondwana and southern Euramerica. For the first scenario, we averaged the number of formations found in eastern and western Gondwana and southern Euramerica for a given time period. The second scenario is similar to the first but included formation counts from northern Euramerica and East Asia in place of western Gondwana. The third scenario only averaged formation counts from eastern Gondwana and southern Euramerica.

 Results Supertree

Topological differences resulted among our supertree, the published source trees, and Marjanović and Laurin's (2019) tree (Fig. 2). In our tree, a polyphyletic “Megalichthyiformes” is the basal-most tetrapodomorph group instead of Rhizodontida (Swartz, 2012 and Zhu et al., 2017). Canowindrids and rhizodontids formed an unexpected sister clade to Eotetrapodiformes. Clack et al.’s (2017) five Tournaisian tetrapod taxa cluster together. Colosteidae is rootward of Crassigyrinus. Caerorhachis is next to Baphetidae. Baphetidae moved crownward compared to previous topologies (likely because of a small character sample size [Marjanović and Laurin, 2019]). Two crownward nodes are unresolved (polytomous). We retained Tungsenia and Kenichthys as the oldest and second oldest tetrapodomorphs. Tristichopteridae, Elpistostegalia, Stegocephalia, Aïstopoda, Whatcheeriidae, Colosteidae, Anthracosauria, Dendrerpetidae, and Baphetidae remain monophyletic. Aïstopoda (Lethiscus and Coloraderpeton) fell rootward to Tetrapoda as reported in Pardo et al. (2017). The average nRF distances quantify differences in topology (Supplementary Table 4). On average, there are 39.7% different or missing bipartitions in the source trees compared to the supertree.

 Phylogeography

We found overwhelming support for a variable rates model of geographic dispersal in early tetrapodomorphs (BF = 632.3; Fig. 3). The estimated rates across the three replicate runs are consistent (out of 122 branches, only three had a median rate scalar with an absolute value difference among the three runs greater than 3). All rate shifts that were two standard deviations greater than the average dispersal rate were reconstructed dispersal events moving from East Asia to southern Euramerica, from eastern Gondwana to southern Euramerica, or southern Euramerica to eastern Gondwana. The fastest estimated dispersal rate occurs along the branch leading to Eotetrapodiformes, moving from eastern Gondwana to southern Euramerica (14.34° × the average rate). As Long et al. (2018) suggest, we find evidence for an East Asian origin for Tetrapodomorpha but with moderate uncertainty (average estimate ± standard deviation of posterior distribution; longitudeavg = 81.5° ± 10.1°, latitudeavg = −6.4° ± 8.5°). We also reconstruct an origin for “Megalichthyiformes” that borderlines East Asia and eastern Gondwana (longitudeavg = 107.2° ± 14.1°, latitudeavg = −22.6° ± 8.7°), along with an eastern Gondwana origin for the clade uniting “Canowindridae” and Rhizodontida (longitudeavg = 137.1° ± 8.2°, latitudeavg = −32.0° ± 4.7°). We recover a southern Euramerican origin for Eotetrapodiformes, consistent with previous studies (longitudeavg = −12.5° ± 7.0°, latitudeavg = −19.4° ± 6.4°). A southern Euramerican origin was also found for Tristichopteridae (longitudeavg = −12.7° ± 6.9°, latitudeavg = −19.7° ± 6.3°) and Elpistostegalia (longitudeavg = −12.3° ± 5.5°, latitudeavg = −13.5° ± 5.3°). As expected in a phylogenetic comparative analysis, uncertainty in estimated node states increases toward the root. However, despite the level of uncertainty within a single run, only three nodes have mean ancestral state values that are greater than an absolute value of 5° among the three replicateruns.

We find good evidence that geographic sampling bias influences dispersal rate estimates, regardless of the route used (95% CI: western Gondwana route U = [800,928]; northern Euramerica + East Asia route U = [832,946]; direct route U = [729,889]; no scenario includes the null U = 499.5; Fig. 4 and Supplementary Figs. 12 and 13). A U-statistic considerably higher than 499.5 suggests that branches with high dispersal rates have lower regional-level formation counts, on average, than branches with low rates. One can also interpret the null U-statistic of 499.5 as a 50% probability that a random branch with a low dispersal rate will rank higher in its regional-level formation count than a random branch with a high dispersal rate. With bootstrapping, we are 95% confident that the probability of a random branch with a low dispersal rate having a higher regional-level formation count than a random branch with a high rate is 72.97–88.99% for the more conservative ‘direct route’ scenario. Under the more liberal ‘northern Euramerica + East Asia route’ scenario, the probabilities are 83.28–94.69%. In sum, branches with high dispersal rates (two standard deviations greater than average) have a smaller number of recorded formations, on average, along their reconstructed dispersal path.

Our results cannot be explained by a fossil record that is more complete through time (Pull of the Recent). A regression model relating regional-level formation count to the minimum age of each branch shows only a weak relationship (slope = −0.044, r 2 = 0.1, P < 0.001). However, the total global (stage-level) formation count (which does not account for geographic variation) does show potential bias from Pull of the Recent (slope = −0.3, r 2 = 0.71, P < 0.0001). If dispersal rates are biased by the increase in number of formations globally, we would also expect to see elevated dispersal rates decrease toward the tips, but a regression model relating stretched branch lengths with time is not supported (slope = −0.025, r 2 = 0.006, P = 0.41).

 Discussion

We expected to infer high dispersal rates for closely related taxa that are distributed across the globe. Our results, unadjusted for geographic bias in the fossil record, confirm this notion. However, we also find a compelling statistical association between high dispersal rates and a low number of formations along dispersal paths—a patchy fossil record is driving inferences of high dispersal rates. Although we did not test for a correlation between dispersal rate and previously used proxies, such as valid taxon count and stage-level formation count, these proxies do not offer clear predictions for explaining dispersal rate variation. High dispersal rate variation is inferred when closely related taxa are geographically separate. For example, valid taxon count cannot explain geographic rate variation because spatial information is lacking in this bias proxy and because sister taxa are likely to have similar counts (these data are phylogenetically structured). Stage-level formation counts will also not explain dispersal rate variation, particularly if high rate variation exists within the same geological stage. Assuming geological formations are evenly exposed and sampled worldwide, low stage-level formation counts should yield geographically variable fossil species and, therefore, drive high dispersal rate variation. However, formations are not evenly exposed or recorded in geological/paleontological databases, including the PBDB. Our formation count table demonstrates this bias (Table 1). Without geographic context, stage-level formation count cannot distinguish between global and local regions. For example, the geological stages that have the highest recorded number of formations are restricted to southern Euramerica where the majority of eotetrapodiform taxa have been discovered. The association between high formation counts in specific regions and high paleobiodiversity in those regions is likely not a coincidence and has a clear impact on how we interpret dispersal history. The earliest tetrapodomorphs are known from China and Australia at geological stages where relatively few formations are recorded outside of East Asia and eastern Gondwana. The basal-most ancestral state estimates reconstruct paleolocations in East Asia (unsurprisingly). This inference (hypothesis) is predicated on the lack of geological formations recorded outside of East Asia during this time period. In addition, the majority of more crownward taxa and their reconstructed ancestral states are located in North America and Europe at geological stages in which relatively fewer formations are known elsewhere. This bias may heavily influence any conclusions made on the location and habitat of the tetrapodomorph water–land transition. Recently discovered taxa could help mitigate this problem by increasing the power of taxon sampling (Heath et al., 2014), such as Tutusius and Umzantsia from South Africa (Gess and Ahlberg, 2018). However, the current lack of cladistic coding for these taxa excludes them from phylogeny-based analyses. The taxonomic resolution of globally-occurring species, like Eusthenodon and Spodichthys, also impacts current models of species dispersal history because of their relatively uniform distribution (Long et al., 2018). Eusthenodon and Spodichthys represent possible cases where taxonomic resolution is too coarse for phylogeographic analyses. Including these species inhibited our MCMC algorithms from reaching convergence. Widely distributed cosmopolitan species that lack intermediate geographic occurrences increase the uncertainty of parameter estimates within phylogeographic models, as is the case here for these two species.

Phylogenetic studies on macroevolution also often fail to incorporate data from the fossil record itself, such as trace fossil occurrences. Non-anatomical data often contribute to our understanding of taxonomic originations, including chiridian (or digit-possessing) tetrapodomorphs for which trace fossil evidence exists about 10 million years before the first elpistostegalian body fossils (Niedźwiedzki et al., 2010). The inclusion of additional data from trace fossils could radically alter our current models of species dispersal history. Finally, it is important to note that the sampling bias proxies are also constrained by database curation biases. Phylogenetic studies on macroevolutionary trends now regularly leverage public databases, such as the PBDB, which allows larger and broader studies. It is unclear how patchy entries, on taxonomic occurrences and geological formations, for example, interact with other biases inherent in the fossil record. Caution is therefore warranted when these databases are mined, as is the case here.

 Conclusions

Phylogenetic studies on macroevolution have not previously incorporated geographic context, which could influence a wide variety of analyses. We demonstrate here that phylogeographic methods are influenced by geographic sampling variability. We develop a simple sampling bias proxy that incorporates geographic information and show that it explains variation in estimated dispersal rates. The majority of elevated dispersal rates are associated with large-scale movements between major landmasses that have very few, if any, relevant geological formations in between. Our analysis is also unlikely to be influenced by “Pull of the Recent”-like effects. Although not the first supertree for early tetrapodomorphs (Ruta et al., 2003), this study presents the first (to our knowledge) with branch lengths, making it useable for phylogenetic comparative analyses. The new supertree comprises many of the major clades previously inferred, but also recovers new ones that will be subject to scrutiny in future studies (discussed further in the Supplementary Material). This supertree should be useful to researchers who aim to use phylogenetic comparative methods to test hypotheses on the evolution of early tetrapodomorphs. In sum, our study estimates ancestral geographical reconstructions consistent with previously hypothesized dispersal patterns in early tetrapodomorphs. We also find that rates of dispersal are strongly influenced by geographic sampling bias. We suggest that researchers incorporate this proxy in phylogeny-based macroevolutionary studies that could be influenced by spatial distribution of the fossil record.

 Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

 Acknowledgements

We thank the MSU Macroevolution Lab, Jack Wilson, and Matt Lavin for helpful discussions, as well as David Marjanović, John Long, and Michel Laurin for their helpful reviews. We thank Ronald and Dee Blakey at Colorado Plateau Geostystems Inc. for granting permission to use the Devonian paleogeographic map and Alida Bailleul and Michel Laurin for helping with the French translations. We also thank Nathalie Bardet, Eric Buffetaut, Annelise Folie, Emmanuel Gheerbrant, Alexandra Houssaye, and Michel Laurin for the invitation to contribute to this special issue celebrating the life and accomplishments of Jean-Claude Rage.

 Supplementary data

Ahlberg P.E. Postcranial stem tetrapod remains from the Devonian of Scat Craig, Morayshire, Scotland Zool. J. Linn. Soc. 122 1998 99–141 10.1006/zjls.1997.0115 Alroy J. Marshall C.R. Bambach R.K. Bezusko K. Foote M. Fürsich F.T. Hansen T.A. Holland S.M. Ivany L.C. Jablonski D. Jacobs D.K. Jones D.C. Kosnik M.A. Lidgard S. Low S. Miller A.I. Novack-Gottshall P.M. Olszewski T.D. Patzkowsky M.E. Raup D.M. Roy K. Sepkoski J.J. Sommers M.G. Wagner P.J. Webber A. Effects of sampling standardization on estimates of Phanerozoic marine diversification Proc. Natl. Acad. Sci. USA 98 2001 6261–6266 10.1073/pnas.111144698 Benson R.B.J. Butler R.J. Uncovering the diversification history of marine tetrapods: ecology influences the effect of geological sampling biases Geol. Soc. Lond. Spec. Publ. 358 2011 191–208 10.1144/SP358.13 Benson R.B.J. Upchurch P. Diversity trends in the establishment of terrestrial vertebrate ecosystems: Interactions between spatial and temporal sampling biases Geology 41 2013 43–46 10.1130/G33543.1 Benson R.B.J. Butler R.J. Lindgren J. Smith A.S. Mesozoic marine tetrapod diversity: mass extinctions and temporal heterogeneity in geological megabiases affecting vertebrates Proc. R. Soc. B 277 2010 829–834 10.1098/rspb.2009.1845 Benton M.J. Donoghue P.C.J. Asher R.J. Friedman M. Near T.J. Vinther J. Constraints on the timescale of animal evolutionary history Palaeontol. Electron. 18 2015 1–106 10.26879/424 Benton M.J. Ruta M. Dunhill A.M. Sakamoto M. The first half of tetrapod evolution, sampling proxies, and fossil record quality Palaeogeogr. Palaeoclimatol. Palaeoecol. 372 2013 18–41 10.1016/j.palaeo.2012.09.005 Budd G.E. Mann R.P. History is written by the victors: the effect of the push of the past on the fossil record Evolution 72 2018 2276–2291 10.1111/evo.13593 Clack J.A. Bennett C.E. Carpenter D.K. Davies S.J. Fraser N.C. Kearsey T.I. Marshall J.E.A. Millward D. Otoo B.K.A. Reeves E.J. Ross A.J. Ruta M. Smithson K.Z. Smithson T.R. Walsh S.A. Phylogenetic and environmental context of a Tournaisian tetrapod fauna Nat. Ecol. Evol. 1 2017 0002 10.1038/s41559-016-0002 Close R.A. Benson R.B.J. Alroy J. Behrensmeyer A.K. Benito J. Carrano M.T. Cleary T.J. Dunne E.M. Mannion P.D. Uhen M.D. Butler R.J. Diversity dynamics of Phanerozoic terrestrial tetrapods at the local-community scale Nat. Ecol. and Evol. 3 2019 590–597 10.1038/s41559-019-0811-8 Coates M.I. Friedman M. Litoptychus bryanti and characteristics of stem tetrapod neurocrania Elliot D.K. Maisey J.G. Yu X. Miao D. Morphology, Phylogeny and Paleobiogeography of Fossil Fishes 2010 Verlag Dr. Friedrich Pfeil München, Germany 389–416 Cope E.D. Synopsis of the Extinct Batrachia of North America Proc. Acad. Nat. Sci. Philadelphia 20 1868 208–221 Criscuolo A. Gascuel O. Fast NJ-like algorithms to deal with incomplete distance matrices BMC Bioinform. 9 2008 166 10.1186/1471-2105-9-166 Criscuolo A. Berry V. Douzery E.J.P. Gascuel O. SDM: A fast distance-based approach for (super)tree building in phylogenomics Syst. Biol. 55 2006 740–755 10.1080/10635150600969872 Didier G. Fau M. Laurin M. Likelihood of tree topologies with fossils and diversification rate estimation Syst. Biol. 66 2017 964–987 10.1093/sysbio/syx045 Didier G. Laurin M. Exact distribution of divergence times from fossil ages and topologies bioRxiv 2018 490003 10.1101/490003 Didier G. Royer-Carenzi M. Laurin M. The reconstructed evolutionary process with the fossil record J. Theor. Biol. 315 2012 26–37 10.1016/j.jtbi.2012.08.046 Dunne E.M. Close R.A. Button D.J. Brocklehurst N. Cashmore D.D. Lloyd G.T. Butler R.J. Diversity change during the rise of tetrapods and the impact of the Carboniferous rainforest collapse Proc. R. Soc. B 285 2018 10.1098/rspb.2017.2730 20172730 Dunhill A.M. Benton M.J. Newell A.J. Twitchett R.J. Completeness of the fossil record and the validity of sampling proxies: a case study from the Triassic of England and Wales J. Geol. Soc. 170 2013 291–300 10.1144/jgs2012-025 Dunhill A.M. Benton M.J. Twitchett R.J. Newell A.J. Testing the fossil record: Sampling proxies and scaling in the British Triassic–Jurassic Palaeogeogr. Palaeoclimatol. Palaeoecol. 404 2014 1–11 10.1016/j.palaeo.2014.03.026 Dunhill A.M. Hannisdal B. Benton M.J. Disentangling rock record bias and common-cause from redundancy in the British fossil record Nat. Commun. 5 2014 4818 10.1038/ncomms5818 Dowding E.M. Ebach M.C. Evaluating Devonian bioregionalization: quantifying biogeographic areas Paleobiology 2019 10.1017/pab.2019.30 Foote M. Origination and Extinction through the Phanerozoic: A New Approach J. Geol. 111 2003 125–148 10.1086/345841 Friedman M. Coates M.I. Anderson P. First discovery of a primitive coelacanth fin fills a major gap in the evolution of lobed fins and limbs Evol. Dev. 9 2007 329–337 10.1111/j.1525-142X.2007.00169.x Gascuel O. Concerning the NJ algorithm and its unweighted version, UNJ Roberts F. Rzhetsky A. Mathematical Hierarchies and Biology 1997 American Mathematical Soc Providence, RI 149–170 Gauthier J. Cannatella D. de Queiroz K. Kluge A.G. Rowe T. Tetrapod phylogeny Fernholm B. Bremer K. Jörnvall H. The Hierarchy of Life 1989 Elsevier Science Publishers B. V. (Biomedical Division) Amsterdam, The Netherlands Gavryushkina A. Welch D. Stadler T. Drummond A.J. Bayesian Inference of Sampled Ancestor Trees for Epidemiology and Fossil Calibration PLoS Comput. Biol. 10 2014 e1003919 10.1371/journal.pcbi.1003919 Gess R. Ahlberg P.E. A tetrapod fauna from within the Devonian Antarctic Circle Science 360 2018 1120–1124 10.1126/science.aaq1645 Heath T.A. Huelsenbeck J.P. Stadler T. The fossilized birth–death process for coherent calibration of divergence-time estimates Proc. Natl. Acad. Sci. USA 111 2014 E2957–E2966 10.1073/pnas.1319091111 Jablonski D. Roy K. Valentine J.W. Price R.M. Anderson P.S. The Impact of the Pull of the Recent on the History of Marine Diversity Science 300 2003 1133–1135 10.1126/science.1083246 Koch C.F. Bias in the Published Fossil Record Paleobiology 4 1978 367–372 Lakner C. van der Mark P. Huelsenbeck J.P. Larget B. Ronquist F. Efficiency of Markov chain Monte Carlo tree proposals in Bayesian phylogenetics Syst. Biol. 57 2008 86–103 10.1080/10635150801886156 Laurin M. The importance of global parsimony and historical bias in understanding tetrapod evolution. Part I. Systematics, middle ear evolution and jaw suspension Ann. Sci. Nat. Zool. 19 1998 1–42 10.1016/S0003-4339(98)80132-9 Lloyd G.T. A refined modelling approach to assess the influence of sampling on palaeobiodiversity curves: new support for declining Cretaceous dinosaur richness Biol. Lett. 8 2012 123–126 10.1098/rsbl.2011.0210 Long J.A. Clement A.M. Choo B. New insights into the origins and radiation of the mid-Palaeozoic Gondwanan stem tetrapods Earth Environ. Sci. Trans. R. Soc. Edinburgh 2018 1–17 10.1017/S1755691018000750 Marjanović D. Laurin M. Phylogeny of Paleozoic limbed vertebrates reassessed through revision and expansion of the largest published relevant data matrix PeerJ 6 2019 e5565 10.7717/peerj.5565 Marshall J.E.A. Reeves E.J. Bennett C.E. Davies S.J. Kearsey T.I. Millward D. Smithson T.R. Browne M.A.E. Reinterpreting the age of the uppermost “Old Red Sandstone” and Early Carboniferous in Scotland Earth Environ. Sci. Trans. R. Soc. 109 2019 265–278 10.1017/S1755691018000968 Niedźwiedzki G. Szrek P. Narkiewicz K. Narkiewicz M. Ahlberg P.E. Tetrapod trackways from the early Middle Devonian period of Poland Nature 463 2010 43–48 10.1038/nature08623 O’Donovan C. Meade A. Venditti C. Dinosaurs reveal the geographical signature of an evolutionary radiation Nat. Ecol. Evol. 2 2018 452 10.1038/s41559-017-0454-6 Pardo J.D. Szostakiwskyj M. Ahlberg P.E. Anderson J.S. Hidden morphological diversity among early tetrapods Nature 546 2017 642–645 10.1038/nature22966 R Core Team R: A Language and Environment for Statistical Computing 2018 R Foundation for Statistical Computing Vienna, Austria Rambaut A. Drummond A.J. Xie D. Baele G. Suchard M.A. Posterior Summarization in Bayesian Phylogenetics Using Tracer 1.7 Syst. Biol. 67 2018 901–904 10.1093/sysbio/syy032 Raup D.M. Boyajian G.E. Patterns of Generic Extinction in the Fossil Record Paleobiology 14 1988 109–125 Revell L.J. phytools: an R package for phylogenetic comparative biology (and other things) Methods Ecol. Evol. 3 2012 217–223 10.1111/j.2041-210X.2011.00169.x Robinson D.F. Foulds L.R. Comparison of phylogenetic trees Math. Biosci. 53 1981 131–147 10.1016/0025-5564(81)90043-2 Ronquist F. Klopfstein S. Vilhelmsen L. Schulmeister S. Murray D.L. Rasnitsyn A.P. A total-evidence approach to dating with fossils, applied to the early radiation of the Hymenoptera Syst. Biol. 61 2012 973–999 10.1093/sysbio/sys058 Ronquist F. Teslenko M. van der Mark P. Ayres D. Darling A. Höhna S. Larget B. Liu L. Suchard M.A. Huelsenbeck J.P. MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space Syst. Biol. 61 2012 539–542 Rothkugel S. Varela S. paleoMap: An R-package for getting and using paleontological maps 2015 Ruta M. Jeffery J.E. Coates M.I. A supertree of early tetrapods Proc. R. Soc. B 270 2003 2507–2516 10.1098/rspb.2003.2524 Sakamoto M. Benton M.J. Venditti C. Dinosaurs in decline tens of millions of years before their final extinction Proc. Natl. Acad. Sci. USA 113 2016 5036–5040 10.1073/pnas.1521478113 Sakamoto M. Venditti C. Benton M.J. ‘Residual diversity estimates’ do not correct for sampling bias in palaeodiversity data Methods Ecol. Evol. 8 2016 453–459 10.1111/2041-210X.2666 Schliep K.P. phangorn: Phylogenetic analysis in R Bioinformatics 27 2011 592–593 10.1093/bioinformatics/btq706 Signor P.W. Lipps J.H. Sampling bias, gradual extinction patterns, and catastrophes in the fossil record Geol. Soc. Am. Spec. Publ. 190 1982 291–296 Stadler T. Sampling-through-time in birth-death trees J. Theor. Biol. 267 2010 396–404 10.1016/j.jtbi.2010.09.010 Swartz B. A marine stem-tetrapod from the Devonian of western North America PLoS ONE 7 2012 e33683 10.1371/journal.pone.0033683 Tennant J.P. Mannion P.D. Upchurch P. Environmental drivers of crocodyliform extinction across the Jurassic/Cretaceous transition Proc. R. Soc. B 283 2016 20152840 10.1098/rspb.2015.2840 Tennant J.P. Mannion P.D. Upchurch P. Sea level regulated tetrapod diversity dynamics through the Jurassic/Cretaceous interval Nat. Commun. 7 2016 12737 10.1038/ncomms12737 Venditti C. Meade A. Pagel M. Multiple routes to mammalian diversity Nature 479 2011 393–396 10.1038/nature10516 Xie W. Lewis P.O. Fan Y. Kuo L. Chen M.-H. Improving Marginal Likelihood Estimation for Bayesian Phylogenetic Model Selection Syst. Biol. 60 2011 150–160 10.1093/sysbio/syq085 Zhang C. Stadler T. Klopfstein S. Heath T.A. Ronquist F. Total-Evidence Dating under the Fossilized Birth–Death Process Syst. Biol. 65 2016 228–249 10.1093/sysbio/syv080 Zhu M. Ahlberg P.E. Zhao W.-J. Jia L.-T. A Devonian tetrapod-like fish reveals substantial parallelism in stem tetrapod evolution Nat. Ecol. Evol. 1 2017 1470–1476 10.1038/s41559-017-0293

Example of how the regional-level formation count proxy is calculated. A. Five major geographic regions are highlighted by color in the Devonian map. Red arrows represent a branch-specific dispersal path to species A, beginning in southern Euramerica and ending in eastern Gondwana. The blue arrow represents the dispersal path to species B. The Devonian map is modified and reproduced with permission from © 2016 Colorado Plateau Geosystems Inc. B. The phylogeny of species A and B scaled by time, with equal branch lengths to both species, and colored to represent the rate of dispersal (red is fast, blue is slow). For every branch of the tree, the number of formations is counted for every region and for each geological stage covered by the dispersal pathway. It is then weighted by the number of geological stages and geographic regions covered. Under the western Gondwana route scenario, the branch to species A covers three geographic regions, while the branch to species B only covers one. Assuming both branches cover only one geological stage, the high dispersal rate for species A can be explained by the lack of recorded geological formations in western Gondwana. C. A line plot of the formation counts through time, colored by geographic region according to the Devonian map above, shows temporal and geographic variability.

Exemple de calcul de l’approximation du nombre de formations au niveau régional. A. Cinq régions géographiques majeures sont figurées en couleur sur la carte du Dévonien. Les flèches rouges représentent une voie de dispersion propre à chaque branche vers l’espèce A, commençant en Euramérique méridionale et se terminant dans le Gondwana oriental. La flèche bleue représente le chemin de dispersion vers l’espèce B. La carte du Dévonien est modifiée et reproduite avec la permission de © 2016 Colorado Plateau Geosystems Inc. B. Phylogénie des espèces A et B extrapolée à l’échelle du temps, avec des longueurs de branches égales pour les deux espèces et colorées pour représenter la vitesse de dispersion (rouge = rapide, bleu = lente). Pour chaque branche de l’arbre, le nombre de formations est compté pour chaque région et est couvert par le patron de dispersion pour chaque étage géologique. Il est alors pondéré par le nombre d’étages géologiques et de régions géographiques couvertes. Dans le scénario de la route du Gondwana occidental, la branche pour l’espèce A couvre trois régions géographiques, tandis que pour l’espèce B, seule une est couverte. Dans l’hypothèse où les deux branches ne recouvrent qu’un seul étage géologique, la vitesse rapide de dispersion pour l’espèce A peut être expliquée par le manque de formations géologiques répertoriées dans le Gondwana occidental. C. Un tracé linéaire rassemblant les comptes de formation à travers le temps, coloré par région géographique, selon la carte du Dévonien ci-dessus, montre la variabilité temporelle et géographique.

The time-scaled tetrapodomorph supertree. Taxonomic groups in quotes are not monophyletic. Here, Glyptolepis, a dipnomorph, is the outgroup. We downloaded the silhouettes from phylopic.org: Eucritta and Greererpeton by Dmitry Bogdanov (vectorized by T. Michael Keesey), Eusthenopteron by Steve Coombs (vectorized by T. Michael Keesey), and Gogonasus and Tiktaalik by Nobu Tamura (CC BY-SA 3.0).

A. Arbre de synthèse des tétrapodomorphes à l’échelle temporelle. Les groupes taxonomiques entre guillemets ne sont pas monophylétiques. Ici, Glyptolepis, un dipnomorphe, est le groupe externe. Les silhouettes ont été téléchargées sur phylopic.org : Eucrita et Greererpeton par Dmitry Bogdanov (vectorisé par T. Michael Keesey), Eusthenopteron par Steve Coombs (vectorisé par T. Michael Keesey) et Gogonasus et Tiktaalik par Nobu Tamura (CC BY-SA 3.0).

A. Trimmed tetrapodomorph phylogeny with mapped rates of dispersal. Cooler (bluish) colors represent slower rates and warmer (reddish) colors represent faster rates. B. Non-eotetrapodiform (left in blue) and eotetrapodiform (right in green) trees and taxon paleolocations plotted on a map of the Middle Devonian. Transparent polygons illustrate broad geographic regions of sampled taxa in southern Euramerica, eastern Gondwana, and East Asia. Numbers show the total number of geological formations recorded from each major geographic region (eastern Gondwana and East Asia combined). Colored circles show average paleolocations of major clades estimated by the ‘geo’ model and indicated in the tree above. Red circle: Tetrapodomorpha, orange: “Megalichthyiformes”, yellow: “Canowindridae” + Rhizodontidae, green: Tristichopteridae, and blue: Elpistostegalia. Phylogeny with mapped dispersal rates was produced in BayesTrees (http://www.evolution.rdg.ac.uk/BayesTrees.html). Middle Devonian tree and paleolocation plots were made using the ‘phylo-to-map’ function in the R package, phytools (Revell, 2012). Middle Devonian map was sourced from the R package, paleoMap (Rothkugel and Varela, 2015). Tetrapodomorph silhouettes were sourced from phylopic.org: Eucritta by Dmitry Bogdanov (vectorized by T. Michael Keesey), Osteolepis by Nobu Tamura, and Acanthostega by Mateus Zica (CC BY-SA 3.0).

A. Phylogénie des tétrapodomorphes et évolution des vitesses de dispersion. Les couleurs les plus froides (bleutées) représentent les vitesses les plus faibles et les plus chaudes (rougeâtres) les vitesses les plus rapides. B. Arbres des non-éotétrapodiformes (à gauche, en bleu) et des éotétrapodiformes (à droite en vert) et paléolaocalisations des taxons rassemblées sur une carte du Dévonien moyen. Les polygones transparents illustrent les vastes régions géographiques des taxons échantillonnés en Euramérique méridionale, au Gondwana oriental et dans l’Est asiatique. Les numéros montrent le nombre total de formations géologiques répertoriées dans chaque région géographique majeure (Gondwana oriental et Est asiatique combinés). Les cercles colorés montrent la moyenne des paléolocalisations des clades majeurs estimés par le modèle « geo » et indiqués dans l’arbre. Cercle rouge : Tetrapodomorpha, orange : « Megalychtyiformes », jaune : Canowindridae + Rhizodontidae, vert : Tristichoptéridae, et bleu : Elpistostegalia. La phylogénie, avec les vitesses de dispersion cartographiées a été produite avec BayesTrees (http://www.evolution.rdg.ac.uk/BayesTrees.html). Les tracés d’arbres et de paléolocalisations du Dévonien moyen ont été réalisés à l’aide de la fonction « phylo-to-map » du R package phytools (Revell, 2012). La carte du Dévonien moyen est issue du R package paleoMap (Rothkugel et Varda, 2015). Les silhouettes de tétrapodomorphes sont issues de phylopic.org: Eucrita par Dmitry Bogdanov (vectorisé par T. Michael Kessey), Osteolepis par Nobu Tamura, et Acanthostega par Mateus Zica (CC BY-SA 3.0).

A. Scatter-plot of the average dispersal rates over the regional-level formation counts for each branch of the phylogeny, using the northern Euramerica + East Asia route scenario. Points colored by the dispersal rate being above (red) or below (blue) two standard deviations greater than the average rate across the tree. B. Histogram of the bootstrapped U-statistics. Values outside of the 95% confidence interval are grayed out. The median and null expected U-statistics are indicated by the red and blue dotted lines, respectively. The null expected U-statistic is based on the null hypothesis that 50% of the branches with low dispersal rates will have greater regional-level formation counts than branches with higher rates. Rejecting the null hypothesis suggests that estimated dispersal rates are biased and correlate with regional-level formation count.

A. Nuage de points des vitesses moyennes de dispersion en fonction du nombre de formations à l’échelle régionale pour chaque branche de la phylogénie , selon le scénario de la route Euramérique nord + Est asiatique. Les points sont en rouge ou en bleu selon qu’ils sont au-delà ou en deçà de deux écarts-types de la vitesse moyenne sur l’arbre. B. Histogramme des statistiques U bootstrapées. Les valeurs hors de l’intervalle de confiance de 95 % sont en grisé. Les statistiques U médianes et la valeur nulle attendue sont indiquées par les lignes pointillées rouges et bleues, respectivement. La statistique Unulle attendue est fondée sur l’hypothèse nulle selon laquelle 50 % des branches à faible vitesse de dispersion auront un plus grand nombre de formations au niveau régional que les branches à vitesse plus élevée. Le rejet de l’hypothèse nulle suggère que les vitesses de dispersion estimées sont biaisées et corrélées au nombre de formations à l’échelle régionale.

Regional- and stage-level (total) formation counts through time. Nombre de formations région et étage géologique.

Table 1 Period Epoch Stage/Age End time (Ma) Northern Euramerica Southern Euramerica Western Gondwana Eastern Gondwana East Asia Total Permian Cisuralian Kungurian 272.95 10 44 11 0 0 65 Permian Cisuralian Artinskian 283.5 2 39 9 0 0 50 Permian Cisuralian Sakmarian 290.1 2 44 4 0 0 50 Permian Cisuralian Asselian 295 2 47 3 0 0 52 Pennsylvanian Late Gzhelian 298.9 1 42 0 0 0 43 Pennsylvanian Late Kasimovian 303.7 0 33 0 0 0 33 Pennsylvanian Middle Moscovian 307 0 16 0 0 0 16 Pennsylvanian Early Bashkirian 315.2 0 28 0 0 0 28 Mississippian Late Serpukhovian 323.2 0 16 0 0 0 16 Mississippian Middle Visean 330.9 0 14 0 1 0 15 Mississippian Early Tournaisian 346.7 0 7 0 1 0 8 Devonian Late Famennian 358.9 1 9 1 5 1 17 Devonian Late Frasnian 372.2 1 11 0 3 2 17 Devonian Middle Givetian 382.7 1 8 1 4 1 15 Devonian Middle Eifelian 387.7 1 8 0 5 2 16 Devonian Early Emsian 393.3 2 5 0 7 3 17 Devonian Early Pragian 407.6 1 6 0 4 2 13 Devonian Early Lochkovian 410.8 1 3 0 3 1 8 Silurian Přídolí Přídolí 419.2 0 0 0 0 1 1 Silurian Ludlow Ludfordian 423 0 0 0 0 2 2 Silurian Ludlow Gorstian 425.6 0 0 0 0 2 2